Kinetic Ising model in an oscillating field: Avrami theory for the hysteretic response 
and finite-size scaling for the dynamic phase transition 
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Hysteresis is studied for a two-dimensional, spin-1/2, 
nearest-neighbor, kinetic Ising ferromagnet in a sinusoidally 
oscillating field, using Monte Carlo simulations and analytical 
theory. Attention is focused on large systems and moderately 
strong field amplitudes at a temperature below T c . In this 
parameter regime, the magnetization switches through ran- 
dom nucleation and subsequent growth of many droplets of 
spins aligned with the applied field. Using a time-dependent 
extension of the Kolmogorov-Johnson-Mehl-Avrami (KJMA) 
theory of metastable decay, we analyze the statistical prop- 
erties of the hysteresis-loop area and the correlation between 
the magnetization and the field. This analysis enables us to 
accurately predict the results of extensive Monte Carlo sim- 
ulations. The average loop area exhibits an extremely slow 
approach to an asymptotic, logarithmic dependence on the 
product of the amplitude and the field frequency. This may 
explain the inconsistent exponent estimates reported in pre- 
vious attempts to fit experimental and numerical data for 
the low-frequency behavior of this quantity to a power law. 
At higher frequencies we observe a dynamic phase transi- 
tion. Applying standard finite-size scaling techniques from 
the theory of second-order equilibrium phase transitions to 
this nonequilibrium transition, we obtain estimates for the 
transition frequency and the critical exponents (/3/u « 0.11, 
j/v ~ 1.84 and v ~ 1.1). In addition to their significance for 
the interpretation of recent experiments on switching in fer- 
romagnetic and ferroelectric nanoparticles and thin films, our 
results provide evidence for the relevance of universality and 
finite-size scaling to dynamic phase transitions in spatially 
extended nonstationary systems. 

PACS number(s): 05.40.+j, 75.60.-d, 77.80.Dj, 64.60.Qb 



I. INTRODUCTION 

The term hysteresis comes from the Greek husterein 
(varepeuj) which means "to be behind" 0]. It describes 
the lagging of an effect behind its cause, as when the 
magnetization of a body lags behind periodic changes 
in the applied field. While the magnetization response 
of a ferromagnet in an oscillating field is probably the 
example most familiar to physicists and engineers ^j-ff), 
hysteresis is a quite common phenomenon. For instance, 



it is also seen in ferroelectrics ||-fl5[ , in which the polar- 
ization lags behind a time-varying electric field. Other 
examples of hysteresis include electrochemical adsorbate 
layers that are driven through a phase transition by an 
oscillating electrode potential in a Cyclic Voltammetry 
experiment |P3|,[^] and liquid-crystalline systems driven 
through a phase transition by pressure oscillations [5j| . 
Recently, a new class of superconducting materials in- 
cluding DyN^P^C have shown hysteresis in the re- 
sistivity when subject to an oscillating magnetic field. 
Acoustic hysteresis in crystals |^0| occurs when the ul- 
trasonic absorption coefficient changes due to oscillations 
in the amplitude of an ultrasonic wave. 

In recent years new experimental techniques, such as 
magnetic force microscopy (MFM) ]2l|-|25[] , have been 
developed that permit measurements of the magnetiza- 
tion state and switching behavior of particles as small as 
a few nanometers. Ferromagnetic particles in this size 
range consist of a single domain in equilibrium. To- 
gether with ultrathin films, they are of interest as po- 
tential materials for ultra-high density recording media. 
The dynamics of magnetization reversal in nanoscale sys- 
tems has been modeled by kinetic Ising systems subject 
to sudden field reversal [^6] . These numerical and 
analytical studies give results in qualitative agreement 
with the experiments mentioned above. Recent experi- 
ments on ultrathin ferromagnetic Fe/Au(001) films |3^ | 
and thin p(lxl) Fe films on W(110) |33] have consid- 
ered the frequency dependence of hysteresis loop areas, 
which were interpreted in terms of effective exponents 
consistent with those found for a continuous spin model 
p2| , p4| -|3^] . Similar experiments have been performed on 
ultrathin Co films on Cu(OOT). A study of this system 
by Jiang et al. |37| reported exponents consistent with a 



mean-field treatment of the Ising model, whereas a recent 
study by Suen et al. |3£| finds very small effective expo- 
nents in the low-frequency regime, apparently consistent 
with the theoretical results we report here. 

The above discussion is far from an exhaustive account 
of hysteresis examples, but it does give an idea of the 
diversity of situations in which this nonlinear, nonequi- 
librium phenomenon is important. Systems that exhibit 
hysteresis have in common a nonlinear, irreversible re- 
sponse which lags behind the applied force. Numerous 
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general mathematical theories have been formulated to 
model hysteretic behavior in a variety of systems, includ- 
ing the Preisach model d]3^-|d| and systems of differ- 
ential equations that display discontinuous bifurcations 
p9| , ^0|j42| . Hysteresis in thermodynamic systems is of- 
ten due to the presence of a first-order phase transition, 
which is the source of strong nonlinearity in the system. 
These details of the nonlinear response can, however, be 
quite different in different systems and even in different 
parameter regimes for the same system. The details must 
be carefully considered in order to accurately predict such 
aspects of the hysteretic response as its dependence on 
the frequency, amplitude, and waveform of the oscillating 
force. Here we present a study of hysteresis in a partic- 
ular model system which incorporates both spatial de- 
grees of freedom and thermal fluctuations and which has 
a first-order equilibrium phase transition. The system 
response in the parameter regime studied in this paper is 
self-averaging and may be described by the Kolmogorov- 
Johnson-Mchl-Avrami (KJMA) theory of metastable de- 
cay gg. 

Specifically, we consider hysteresis in a spin- 1/2, 
nearest-neighbor, kinetic Ising ferromagnet on a two- 
dimensional square lattice with periodic boundary condi- 
tions, which is subject to a sinusoidally oscillating field. 
For convenience, and because many of the experimen- 
tal measurements of hysteresis involve magnetic systems, 
we use the customary magnetic language in which the 
order parameter is the dimensionless magnetization per 
site, m{t) £ [—1, +1], and the force is the magnetic field 
H(t). However, we expect our results also to apply to 
hysteresis in other areas of science. For example, in di- 
electrics m(t) and H{t) can be re-interpreted as polariza- 
tion and electric field, in adsorption problems as coverage 
9(t) = [m{i) + 1] /2 |4(| and (electro)chemical potential 
or (osmotic) pressure, etc. 

Below its critical temperature and in zero field, this 
model has two degenerate ordered phases corresponding 
to a majority of the spins in the positive or the nega- 
tive direction. A weak applied field breaks the degener- 
acy, and the phase with the spins aligned (anti-aligned) 
with the field is stable (metastable). If the field varies 
periodically in time, the system is driven back and forth 
across a first-order phase transition at H — 0, and the two 
phases alternate between being momentarily stable and 
metastable. As a result, m(t) lags behind H(t), and hys- 
teresis occurs. In the regime of large system size, mod- 
erately strong field, and temperature well below T c con- 
sidered here, the system switches smoothly and almost 
detcrministically between the two magnetized phases. 

The metastable phase in Ising models subject to a sud- 
den field reversal from H to — H decays by different mech- 
anisms, depending on the magnitude of H, the system 
size L, and the temperature T. Two distinct regimes are 
separated by a crossover field called the dynamic spin- 
odal, i?DSP ~ (lnL)^ 1 ^ d ^ 1 \ where d is the spatial di- 
mensionality |l7|,[48| • Detailed discussions of these differ- 
ent decay modes are found in Refs. f^,|§|4|]. At suf- 



ficiently low T that the single-phase correlation lengths 
are microscopic, the different decay regimes can be dis- 
tinguished by the interplay among four length scales: the 
lattice spacing a, the system size L, the radius of a critical 
droplet R c oc l/|iJ|, and the average distance a supercrit- 
ical droplet interface propagates before encountering an- 
other droplet Rq cx exp{S (T)/ [(d + l)!^" 1 ] }. The 
physical significance of 3q(T) is explained in Sec. |n]. In 
this paper, we specifically consider decay in the multi- 
droplet (MD) regime [H DSP {T,L) < \H\ < H MFSP (T) 
where R c (Hmfsp) ~ 0.5]. In terms of the characteristic 
lengths, the MD regime is defined by 



a < R c < R < L 



(1.1) 



Here, the decay of the metastable phase proceeds by ran- 
dom homogeneous nucleation of many critical droplets of 
the stable phase, which then grow and coalesce. The MD 
regime is distinct from the strong-field regime, \H\ > 
Hmfsp(T). It is also distinct from the single-droplet 
(SD) regime (\H\ < H BSP (T, L)), where a < R c < L < 
R . In the SD regime, the decay of the metastable phase 
proceeds by random homogeneous nucleation of a single 
critical droplet of the stable phase. Hysteresis in that 
regime is described in detail in Ref. [p0|. The present 
study, as well as our previous work [ |51| , |52| |, shows that 
the response to an oscillating field is significantly differ- 
ent in the MD and SD regimes. 

Theoretical and computational studies of hysteresis 
have been performed for several models, using a variety 
of methods |Q . These include various studies of models 
with a single degree of freedom, equivalent to mean-field 
treatments of the Ising model |54|-|56||, Monte Carlo (MC) 
simulations of the spin- 1/2 Ising model [ ^5| , p6|57| |66|, and 
several O(N) type models |L2],|| HIlF Thcsc stu dies 
were performed with variations in the details of the sim- 
ulations and in the model parameters. Most of them 
indicate that the average hysteresis-loop area appears to 
display power-law dependences on the frequency and am- 
plitude of H(t). However, there is no universal agreement 
on the values of the exponents, either experimentally or 
theoretically. For the Ising model it has been pointed out 
that nucleation effects would lead to an asymptotically 
logarithmic frequency dependence |ll| , p9| , p5| , p8| . A mean- 
field model exhibits a dynamic phase transition in which 
the mean period-averaged magnetization changes from a 
nonzero to zero mean value p5[ . Such a dynamic phase 
transition has been observed in MC simulations of a ki- 
netic Ising model as well [^7|- |62| , |6~l| [36l|6S|l . A fundamen- 
tally different example of criticality in a hysteretic system 
is the zero-temperature, random-field Ising model, which 
exhibits critical behavior in the hysteresis loop as a func- 
tion of disorder Q . 

The work 
presented in this paper and in Refs. |50-5^j6^] differs 
from most past theoretical and numerical studies of hys- 
teresis in two important ways. First, mean-field models 
do not take into account thermal noise and spatial vari- 
ations in the order parameter, thus ignoring fluctuations 
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which may be important in real materials. Second, most 
previous investigations of hysteresis in Ising models have 
considered the frequency and amplitude dependence of 
quantities such as the loop area and the period-averaged 
magnetization without considering the manner in which 
the metastable phase decays. In this paper, the long- 
time behavior of the hysteretic response is analyzed by 
studying the power spectral densities of the magnetiza- 
tion time series as well as the statistical properties of 
the period-averaged magnetization Q, the loop area A, 
and the correlation B. These quantities arc defined as 
follows. 



Q = — <p m(t) dt 
2ir J 



A 



B 



2tt 



m(H) dH 



m{t) H{t) dt , 



(1.2) 
(1.3) 
(1.4) 



where the initial time t of the period is defined such that 
H(t ) = H(t + 2-Kjuj) = 0. 

Due to the multi-droplet decay mechanism, the average 
hysteresis-loop area exhibits an extremely slow crossover 
to a logarithmic decay with frequency and amplitude in 
the asymptotic low-frequency limit |52j| . This asymptotic 
behavior of the loop area in the MD regime is qualita- 
tively similar to that for the SD regime for two dimen- 
sions [^0 52 1 . However, the calculation is somewhat more 
involved, and the quantitative behavior of the loop area 
is different for the two regimes. 

Beside our results on the low-frequency loop areas, our 
most significant finding is detailed evidence of finite-size 
scaling at a dynamic phase transition (DPT) in the MD 
regime. Here we provide a full account of these results, 
which we briefly reported in Ref. fl69fl . This transition 
can be intuitively understood as a competition between 
two time scales: the period of the external field, 2tt/u>, 
and the average lifetime of the metastable phase, (t(H)), 
defined as the first-passage time to a magnetization of 
zero following an instantaneous field reversal from H to 
-H. If 2tt/lu < (t(H )) [H is the amplitude of H(t)\ 
the magnetization cannot fully switch sign within a sin- 
gle period, and \Q\ > 0. We shall refer to this situation 
as the ordered dynamic phase. If 2n/uo (t(Ho)) the 
magnetization follows the field, and Q « 0. This is the 
disordered dynamic phase. Between these limits there 
is a critical frequency at which appears to become 

singular in the infinite-system limit. We emphasize that 
the DPT is a noncquilibrium phase transition and that 
the probability distribution of the system magnetization 
which characterizes the two phases never relaxes into a 
stationary state. However, the "filtered" time series of Q 
for successive field periods is a stationary stochastic pro- 
cess. To avoid confusion we establish the following ter- 
minology. By the term "dynamic phase" we mean one of 
the qualitatively different system responses separated by 
the DPT. In contrast, the term "phase" by itself refers in 



the conventional way to a uniform thermodynamic phase. 

The rest of this paper is organized as follows. Details 
of the model and a brief review of relevant aspects of the 
Avrami theory of metastable decay are given in Sec [fi|. 
Time-series data for the magnetization and the period- 
averaged magnetization Q are discussed in Sec. |l|. In 
Sec. IV we discuss the power spectral densities obtained 
from the time-series data. In Sec. [v| we obtain an analyt- 
ical result for the time-dependent system magnetization 
during a single period of the field, based on the droplet 
theory of nucleation and a time-dependent extension of 
the Avrami theory. Section [v| contains an analysis of the 
hysteresis-loop area A and the correlation B. This sec- 
tion includes MC data for the probability distributions 
and averages of A and B along with theore tical predic- 
tions based on the results in Sec. In Sec. |VII| we con- 



sider the period-averaged magnetization Q, which is the 
order parameter for the dynamic phase transition. De- 
tails of the finite-size scaling analysis for this transition 
are given in this section as well. A summary, discussion, 



and topics for future study are presented in Sec. VIII 



II. MODEL 

The model used in this study is a kinetic, nearest- 
neighbor Ising ferromagnet on a hypercubic lattice with 
periodic boundary conditions. The Hamiltonian is given 
by 



H 



-jY^s l s ] -H{t)Y J 



(2.1) 



(ij) 



where H (t) = — Ho sm(uit) , Si = ±l is the state of the ith 
spin, Yluj) runs over all nearest-neighbor pairs, and ^ 
runs over all N — L d lattice sites. The magnetization per 
site is 



m{t) 



1 N 



(2.2) 



The dynamic used is the Glauber J7l| single-spin-flip 
Monte Carlo algorithm with updates at randomly chosen 
sites. The time unit is one Monte Carlo step per spin 
(MCSS). The system is put in contact with a heat bath 
at temperature T, and each attempted spin flip from Si 
to —Si is accepted with probability 72 



W{ Sl 



exp(-pAEj) 
1 + exp(-f3AE t ) 



(2.3) 



Here AEi is the change in the energy of the system that 
would result if the spin flip were accepted, and j3= l/k^T 
where fee is Boltzmann's constant. It has been shown in 
the weak-coupling limit that the stochastic Glauber dy- 
namic can be derived from a quantum-mechanical Hamil- 
tonian in contact with a thermal heat bath modeled as 
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a collection of quasi-free Fermi fields in thermal equilib- 
rium [fr3f . 

The average number of droplets of the stable phase 
that are formed per unit time and volume is given by the 
field and temperature dependent nucleation rate per unit 
volume, 



related. In the simple field-reversal case this leads to the 
well known "Avrami's Law," 



I(H(t),T)*B(T)\H(t)\ K exp 



SoCO 



\H(t)\ 



d-i 



(2.4) 



The notation follows that of Ref. j26|, where B(T) is a 
non-universal temperature dependent prefactor, and K 
and Sq(T) are known from field theory |j74|,|?5), MC sim- 
ulations Eq] and numerical transfer-matrix calculations 
|^6|,0; their values are listed in Table |. The quantity 
So(T) is the field- independent part of the free-energy cost 
of a critical droplet, divided by fc^T. The field, H(t), is 
the only quantity through which I (H(t), T) depends on 
time in this adiabatic approximation. 

Several other quantities, whose values do not depend 
on the frequency of the field, are required as input for the 
theoretical calculations in the following sections and are 
listed in Table |. They are determined through what we 
refer to as "field-reversal simulations." In these simula- 
tions the system initially has all spins up, i.e. positive. It 
is then subjected to a static field of magnitude Hq with 
a sign opposite the system magnetization. This instan- 
taneous field quench prepares the system in a metastable 
state, and the decay of this metastable phase proceeds 
by the MD mechanism outlined in the introduction. 

In the MD parameter regime the average size of a criti- 
cal droplet R Cl and the average distance between droplets 
Ro, are both much smaller than the system size. There- 
fore, many droplets nucleate and grow to drive the system 
into the stable phase, resulting in an almost determinis- 
tic decay process. This can be understood by imagining 
the system subdivided into cells of linear size Rq 7 which 
each contain a single droplet. Each of these subsystems 
will appear to be in the SD regime, and the time taken 
to nucleate a critical droplet is stochastic with an ex- 
ponential probability distribution. However, as a con- 
sequence of the Central Limit Theorem the probability 
density of the lifetime for the entire system asymptoti- 
cally approaches a Gaussian as the system size increases. 
The decay process thus becomes increasingly determinis- 
tic, with a lifetime distribution whose variance decreases 
as L~ d with increasing system size ]26| . 

In addition to the "self-averaging" process described 
above, two other concepts are needed to understand the 
MD decay process: the droplet interface velocity and the 
overlapping of growing supercritical droplets. Our de- 
tailed treatment of these effects is given in Sec. [v], where 
the theory of MD decay in static field is generalized to 
time-dependent fields. Accounting for all of these effects, 
an expression may be obtained for the time-dependent 
magnetization of a system in a field-reversal experiment. 
In the KJMA approximation p^-p5|], it is assumed that 
the positions and sizes of the growing droplets are uncor- 



i(t) = 2 exp [-$(*)] - 1 



2 exp 
2 exp 



m d {v t') d dt 



d+1 



(2.5a) 
(2.5b) 

(2.5c) 



where vo, the interface velocity for a growing droplet 
of the stable phase, is assumed to be constant, Q d is 
a proportionality constant such that the volume V of a 
droplet of radius R is V = £l d R d , and the other constants 
have been introduced previously. The integral $(t) in 
Eq. (2.5b) is the "extended volume" p5| , i.e. the total 
volume fraction of droplets of the equilibrium phase at 
time t, uncorrected for overlaps. The assumption that the 
positions and sizes of the droplets are uncorrelated leads 
direc tly t o the exponential relation for m(t) |78|. Solving 



Eq. (2.5c) for the time at which m = gives the average 



lifetime in the MD regime, 



n d vp 



>2(d + l) 



(2.6) 



which, in contrast to the SD regime, is independent of L. 
To describe the hysteretic response in the MD regime in 
Sees. ^ and [v| we employ a time-dependent extension of 
this theory, in which / and vq are both functions of time 
through H(t). 



III. TIME-SERIES DATA 

All numerical simulations reported in this paper are 
performed for d = 2, T = 0.8T C and one of three system 
sizes, L = 64, 90, and 128. A sinusoidal field is applied 
to the system with amplitude, Ho = 0.3J > Hbsp(L), 
chosen such that in field-reversal simulations the system 
is clearly in the MD regime for a field of magnitude Hq 
for all three values of L. This is illustrated in Fig. [y. 
The dynamic spinodal field is approximated by i?DSP ~ 
if 1/2, where H 1 / 2 is the value of H (for given L and T) 
for which the relative standard deviation of the lifetime, 
r = (7 T / (t) is 1/2 (oy is the standard deviation of the 
lifetime). The values of H DSP and (r> for T = 0.8T C 
and L — 64 |7?||8(J are given in Table B. For the larger 
systems, (r) is approximately unchanged while -Wdsp and 
r are smaller than for L = 64. Thus the system is well 
within the MD regime for all three sizes used. 

To obtain the raw time-series data, the system was 
initially prepared with either a random arrangement of 
up and down spins with m(t = 0) « 0, or with a uni- 
form arrangement with all spins up. Then the sinusoidal 
field was applied and changed every attempted spin flip, 
allowing for a smooth variation of the field. The time 
series did not appear to depend on the initial conditions 
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after a few periods. For each system size, the simula- 
tions were performed with several values of the driving 
frequency uj. For each frequency, we recorded the time- 
dependent magnetization m(t). Most of the simulations 
at intermediate and high frequencies were recorded for 
approximately 16.9 x 10 6 MCSS. (Simulations for some 
of the lowest frequencies were recorded for approximately 
5.9 x 10 5 MCSS.) Files containing the data for these 
longest runs are about 800 megabytes and required 9 
days (one month) to run for L = 64 (£ = 128) on a single 
66 MHz node of an IBM sp2 computer. Since the hys- 
teresis depends on the competition between the two time 
scales represented by the field period and the metastable 
lifetime, we chose the frequencies of H(t) by specifying 
the ratio 



R = 



(2vr/cj) 

Vim 



(3.1) 



One may think of R as a scaled period and 1/R as a 
scaled frequency. 

We note that (t(H )) is the "shortest of the long time 
scales" that describe the system. For T = 0.8T C and all 
the values of L used here, the time scale for spontaneous 
fluctuations between the phases in the absence of an ap- 
plied field, (t(0)), is essentially infinite. Even when the 
field has its maximum strength Hq, the nucleation of the 
critical droplets necessary to leave the metastable phase 
is always driven by the thermal fluctuations. Driving the 
system from the metastable to the stable phase therefore 
truly depends on the joint action of the random ther- 
mal noise and the deterministic oscillating field. Figure 
U shows short initial segments of the magnetization time 
series for four different values of R. The figures for the 
first three values of R are chosen to represent: a sys- 
tem in the ordered dynamic phase (R = 3, (\Q\) > 0), 
near the dynamic transition (R — 3.436 « -Rcr), and in 
the disordered dynamic phase (R = 7, (\Q\) « 0). This 
value for R CI is obtained by finite-size scaling analy sis o f 



the probability density for Q, as described in Sec. [VI] . 
The time-series segment shown for i? = 200 is deep in the 
disordered dynamic phase region and illustrates the be- 
havior of the system for very low frequencies. The stan- 
dard deviation of the average lifetime in the MD regime 
is relatively small compared to the SD regime. If the pe- 
riod of H(t) is sufficiently long, the system has enough 
time to switch phases during a single half-period. This is 
clearly seen at R = 7, for which the system switches dur- 
ing practically every half period. Soon after m(t) reaches 
saturation, the field passes through zero and favors the 
opposite phase. Similar behavior is seen in the time se- 
ries for R = 200, except that the period is so long that the 
system decays to the favored phase well before the field 
reaches its maximum value. Then the magnetization fluc- 
tuates near its equilibrium value m(t) ±1 until the field 
again switches sign and the system once more becomes 
metastable. If the period of H(t) is sufficiently short, 
the system does not have time to switch during a single 
half-period. This can be seen for i? = 3. While the field 



favors the opposite phase, the magnetization changes as 
many critical droplets nucleate and begin to grow. Be- 
fore the magnetization can completely reverse however, 
the field changes sign and the droplets of the now un- 
favored phase shrink and disappear. For R — 3.436, the 
period near the critical value, the period-averaged mag- 
netization slowly "meanders" from positive to negative 
values over several periods of H(t). This "slow switch- 
ing" occurs many times over the entire time series. The 
number of field cycles shown in Fig. |2| is small compared 
to the total number of cycles in an entire time series. 

The "slow switching" seen near the dynamic phase 
transition also occurs for frequencies in the ordered dy- 
namic phase region, where (|Q|) > 0. However, there the 
times between consecutive switches are too long to show 
in plots of m(t) vs t. For this reason, the "filtered time 
series" for Q are shown in Fig. ^, which provide plots 
of Q for consecutive periods in the magnetization time 
series. Even for the Q time series, the number of periods 
shown is small compared to the total number of periods, 
except for i? = 200 which displays the entire time series. 
For the low frequencies, the Q values are concentrated 
near Q = 0, with larger fluctuations for R = 7 than for 
i? = 200. Analysis of the Q data and the dynamic phase 



transition will be detailed in Sec. VII 



IV. POWER SPECTRAL DENSITIES 

A standard method used to characterize a time series 
is to calculate its power spectral density (PSD). Figure 
^ shows the PSDs of the raw data, short segments of 
which are shown in Fig. ^| The technical details on how 
the PSDs were obtained are elaborated in Ref. ]5(J. For 
clarity, the PSDs for different driving frequencies, shown 
in Fig. |](a) with L — 64, have been shifted in the verti- 
cal direction by arbitrary offsets. The same spectra are 
plotted in Fig. 0(b) with no offset. The fourth spectrum 
shown in Fig. |[, labeled "background," corresponds to 
thermal equilibrium fluctuations in a single thermody- 
namic phase. To obtain this spectrum, a simulation was 
performed on a system with the same size, temperature, 
and for the same number of MCSS as the other spectra, 
in a static field of 

To describe the PSD for each frequency, we identify 
three distinct regions: 1) the peaks, 2) the thermal noise 
region, and 3) the low-frequency region. The most promi- 
nent features of the PSDs are the sharp peaks. For R = 3 
and 3.436 ~ i? cr , the first peak in the spectrum corre- 
sponds to lu, the frequency of the external field. The 
second peak corresponds to 2ui, and so on. These odd 
and even harmonic peaks arise because the shape of the 
time series is not purely sinusoidal due to the nonlin- 
ear response of the system. For R~7 (and longer pe- 
riods), only odd harmonics are seen. The extinction of 
the even harmonic peaks occurs because the shape of 
the time series is beginning to resemble a square wave 
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(see Fig. g). The power p n contained in the nth com- 
ponent of the Fourier series for a pure square wave is 
p n = 16[sin(n7r/2)] 4 /(wr) 2 , which decays as n~ 2 and 
vanishes for even n. This enables one to understand the 
reduced second harmonic in the PSD for R = 7, which is 
just barely observable between the first two sharp peaks. 
However, in contrast to our observations in the SD regime 
|pOfl , no dips in the PSD at even n, corresponding to the 
zeros oi p n |81|, were observed for the values of R ana- 
lyzed here. 

Unlike the SD regime |Q, the highest frequencies for 
each of the PSDs do not fall onto the thermal noise back- 
ground. Since the average lifetime in the MD regime is 
much smaller than in the SD regime, the R values shown 
in Fig. |] correspond to much larger frequencies in units 
of MCSS" 1 . Therefore, the time scales characterizing the 
peaks and the thermal noise regions are not as well sep- 
arated as for the SD regime, so the two regions overlap. 
Also, none of the PSDs shown here are for low enough 
frequencies that the metastable phase can decay and the 
system remain in the stable state sufficiently long during 
each half-period to sample the purely thermal fluctua- 
tions which would display exponential time correlations. 

The low-frequency region comprises the portion of 
each spectrum between the first peak and the lowest 
resolved frequency. The PSD in this region exhibits a 
strong dependence on the frequency of the field. Signif- 
icant amounts of power in this portion of a PSD indi- 
cate the presence of slow behavior on time scales larger 
than that of the driving field. Near the transition fre- 
quency, R — 3.436 f» i? cr , the overall slope in the low- 
frequency region is close to —2. This suggests that the 
long-time correlations are exponential with a very large 
correlation time |8^| . The turnover in the corresponding 
Lorentzian PSD is not observable because of the lack of 
low-frequency resolution due to the finite length of the 
time series. For R — 3, the low-frequency region of the 
PSD also suggests a Lorentzian, again with a correlation 
time that is difficult to estimate due to the poor low- 
frequency resolution. For R = 7, the flat low-frequency 
region is that of the PSD of white noise. This is con- 
sistent with the behavior of the Q time series shown in 
Fig. ||(d). The PSDs for other system sizes display a 
qualitatively similar frequency dependence in the sharp 
peak and low-frequency regimes discussed above. How- 
ever, there is a systematic size dependence in the PSDs 
for the thermal noise background, which is smaller for 
larger L. This is easily understood since the variance of 
the magnetization in equilibrium should be proportional 
to L- 2 . 



V. DERIVATION OF M(T) FROM AVRAMI'S 
LAW 

The theoretical predictions for the frequency depen- 
dence of both the hysteresis-loop areas and the cor- 
relation rely on numerically solving an expression for 



m(t) during a single field period. This expression is de- 
rived from a generalization of Avrami's law with time- 
dependent nucleation rate and growth velocity. Using 
these values for m{t) in Sec. VI we explicitly calculate the 
average hysteresis-loop area, (A) = — (<f m{H) dH), and 
the average correlation, (B) = (u>/2tt) (§ m{t) H(t) dt). 
Avrami's law gives the volume fraction of the metastable 
state for systems in which many non-interacting droplets 
nucleate, then grow and coalesce without changing their 
shapes as the system decays to the stable phase. The vol- 
ume fraction of metastable phase is related to the magne- 
tization by cj)(t) « [m(t) + 1] /2. The original KJMA cal- 
culation shows that the volume fraction of the metastable 



phase, 4>{t), is given by Eq. (2.5a). For a homogeneous, 
time -depe ndent nucleation rate, the extended volume in 
Eq. (|2.5h| ) is generalized to 



$(*)= f I(L,T;t n )V(t,t n )dt B , 
Jo 



(5.1) 



where V(t, t n ) is the volume at time t of a droplet which 
was nucleated at time t n , and I(L,T;t n ) is the time- 
dependent nucleation rate given by Eq. (2.4). The vol- 
ume V(t, t n ) is given by 



V(t,t n ) 



t -id 

v(t')dt' 



(5.2) 



The Lifshitz-Allcn-Cahn approximation |8J| |85|] is used 
to specify the interface velocity of a growing droplet as 
v(t) w v\H(t)\. (The proportionality factor v should not 
be confused with the critical exponent discussed in 
Sec. VII.) The effects of the dependence of the propor- 



tionality factor v on the droplet radius, which are re- 
mar kablyminor, are discus sed by Shncidman and collab- 
orators l86l,p7|. Equation (5.1) with d = 2 and A' = 3 in 



Eq. (|2.4|) then gives 



dt n (5.3) 



$(i)=^ I(L,T;t n )Q 2 uH smut'dt' 

pt q v 2 H 2 

= / I(L,T;t n ) — - [cos tut n — cos tot] dt n 

Jo ^ 

B(T)Q, 2 v 2 H?: /"* r l2 . 

= 5 / [cos ujt n — cos ujt\ \Hosmcjt n \ 

^ Jo 



x exp 



'o 

seen 



| Hq sin tot n | 



dt n . 



Numerically integratin g Eq . (5.3) for t 6 [0,n/w] and 
substituting into Eq. (2.5a) enables one to find <j>(t) = 
exp [— $(t)] for a system which starts with c/!>(0) = 1. The 
only quantity which must be determined from the MC 
simulations is B(T), the field-independent part of the 
prefactor in the nucleation rate. In the present paper, 
this single parameter is set by demanding that the aver- 
age simulated loop area for i? = 200 matches the theoret- 
ical prediction (see Sec. VI). 
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The integration of Eq. ( |5.3| ) can only be performed for 
the first half-period of the field cycle, in which m(t) and 
H(t) have opposite signs. When m(t) and H(t) have 
the same sign, the droplets that were formed during the 
previous half-period will shrink. In that case, the mag- 
netization is unable to switch sign during a single pe- 
riod for sufficiently high frequencies (see Fig. ||(a)). The 
KJMA theory cannot be used to find the volume fraction 
of shrinking droplets of the metastable phase. Therefore, 
an approximate calculation of m(t) for t € [ir/u), 2tt/uj] 
must be devised. Assume that the shrinking droplets in 
Fig. ^|(c) are described by merely reversing the nucleation 
and growth process represented in Fig. ||(b). Then, the 
volume fraction of the growing, now stable, background 
for t € [k/o->, 2tt/uj] would be given by </> [(2ir/ui) — t]. In 
addition to the growing stable background, which en- 
croaches upon the droplets of the metastable phase dur- 
ing the second half period, we assume that droplets of 
the stable phase nucleate and grow within the shrinking 
metastable droplets as depicted in Fig. ||(c). Further- 
more, we assume that this nucleation and growth pro- 
cess is also described by Avrami's law, with the shrink- 
ing metastable phase analogous to the metastable back- 
ground filling the entire system for t £ [0, 7r/o;]. The com- 
plete prescription for the volume fraction of the growing 
stable background for t £ [tt/uj, 2ir/ui] is then 



4>'{t) « 



t*J + J [shrinking metastable phase] 



x [growing stable phase inside shrinking droplets] | 

Id 

2tt 



2tt 



t + 



1 



t 









'(«-=)] 







(5.4) 



Thus the theoretical result for the magnetization during 
an entire period of H (t) is 



1 



m(t) 



2<f>'(t) - 1 = 1 - 2(f) i 



< t < tt/uj 

tt/lo < t < 2tt/uj . 

(5.5) 



Note that, in this approximation m{2ir/uj) = m(0) for all 
frequencies. Although it gives very good values for (A) at 
all frequencies and for (B) especially at low frequencies, 
it gives a continuous change in (Q) with to, with no sign 
of a dynamic phase transition. 

It is possible to rewrite Eq. ( |5.3|) for $(£) by making 
the substitution H = Hosin(ujt) and H' = Hosm(ojt n ). 
After tedious, but straightforward algebra this gives 

m) = B ( T )^f H o [" H , 3e s oiT)/H > 



I ( j£ 1 - 12 V I - 



1 1 



+ 



dH' 



The equation is valid for H e [0, Ho], so this expres- 
sion allows a theoretical calculation of the magnetization 
for sufficiently low frequencies such that m(t) completely 
switches sign within a quarter period of the field. Notice 
that oj is no longer part of the integration operation. A 
similar calculation has been used to describe experimen- 
tal data for switching currents in (Nb,Co)-doped BaTiC>3 
ceramics |l3| and hysteresis loops for triglycine sulphate 
(TGS) single crystals (ljjl^] in the presence of an os- 
cillating electric field. Those studies assume heteroge- 
neous nucleation and therefore obtain a frequency de- 
pendence of ui~ d , rather than 

uj-id+i) [ see Eq . Q] f or 

the case of homogeneous nucleation studied here. The 
integral in Eq. ( |5.6D is evaluated once for H e [0, H ], 
from which m(t) for t S [0, 7r/(2w)] at any frequency can 
be calculated. For numerical reasons, Eq. (5.6) enables 



the theoretical prediction of the hysteretic response to 
be more e asily calculated for lower frequencies than by 
using Eq. (5.3) alone. The values of m(t) obtained from 
both Eq. ( |5.3| ), for intermediate to high frequencies, and 
Eq. ( |5.6| ), for low frequencies, are used in calculating the 
theoretical predictions for (^4) and (B) in Sec. VI. 



VI. 



HYSTERESIS-LOOP AREAS AND 
CORRELATION 



In this section we calculate the hysteresis-loop area 
A and the correlation between the field and the system 
magnetization B, defined in Eqs. (1.3) and (1.4) respec- 
tively. These quantities are calculated over each period 
in the entire time series. From the resulting "filtered" 
time series we construct histograms to obtain the prob- 
ability densities of A and B for each separate frequency 
of the field. The hysteresis-loop area represents the en- 
ergy dissipated during a sing le period of the field §-|. 
It is therefore one of the most important physical quanti- 
ties characterizing hysteretic systems, and it is frequently 
measured in experiments. 

Recent experiments on ultrathin ferromagnetic films 
p2| , p3| , as well as numerical simulations of two- 
dimensional Ising models [ p5p^j6^ , |8^ ] , have been inter- 
preted in terms of a low-frequency power law, A oc H^uj b , 
with a range of exponent values reported [ |62p^j8^] . This 
interpretation is not fully consistent with the fluctuation- 
free mean-field result ]54] , p6| ], A = Aq + const [uj 2 (Hq — 
H^p)] 1 / 3 with positive constants A and H sp , which has 
been applied to analyze experiments on ultrathin films 
of Co on Cu(001) |7J. Nor does the single power-law 
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dependence agree with the logarithmic dependence ex- 
pected if thermally activated nucleation JlJ,|2^,^,^8[ is 
the rate-determining process. Here we present in detail 
analytical and numerical results that indicate a resolu- 
tion of this puzzling situation. A brief summary of some 
of the results was given in Ref . j52| . 

The hysteresis loops depend on frequency and have 
qualitatively different shapes above, near, and below the 
dynamic phase transition. Figure [3] shows plots of the 
magnetization vs field in the MD regime for four field fre- 
quencies. The loops shown here are for the same frequen- 
cies as the time series shown in Fig. |[ The loops shown 
in Fig. |^(b) indicate the large fluctuations in the average 
magnetization near the transition at R = 3.436 « R cr . 
As the frequency gets further away from the transition, 
the variation in the loops becomes smaller. The loops for 
a very low frequency, i? = 200 in Fig. ^|(d), show how the 
magnetization reverses sign early in each half-period of 
the field, saturating close to the equilibrium magnetiza- 
tion at a field well below Hq. The low- frequency approx- 
imations for A derived later in this section are applicable 
to these squarish loops, in which the switching field is 
small and the field variation is nearly linear. 

Probability densities for A are shown in Fig. 7| The 
loop area has been divided by the maximum possible loop 
area, AHq- Figure 0(a) shows the probability densities for 
L = 64 and several values of R between R = 2 and R = 200. 
In contrast to the SD regime (Ho), the distributions are 
unimodal for all frequencies. The roughness of the dis- 
tributions at the lowest frequencies is due to the lack 
of statistics because the time series for these frequencies 
contain a smaller total number of periods than the data 
for the intermediate and high frequencies. Figure f7|(b) 
shows the probability densities for a frequency near the 
transition for L = 64, 90, and 128. The position of the 
peak for each of these distributions appears to be inde- 
pendent of system size. However, the width of the distri- 
bution becomes smaller with increasing system size. This 
observation is consistent with the system size dependence 
of the switching-field distribution in field-reversal exper- 
iments ||(|. For further details on the ^-dependen ce o f 
the width near the transition, see the end of Sec. VH. 



for t s becomes 



The average loop area for a specific frequency is a quan- 
tity often displayed in experimental and numerical stud- 
ies of hysteresis. It can also be obtained in the present 
case, and we do so below. The means of the distributions 
in Fig. 0(a) are shown in Fig. along with the average 
loop areas for L = 90 and 128. The solid cu rve in Fig. |]is 
obtained by numerical integration of Eq. (1.3) using the 
values of m(t) obtained in Sec. |v|. 

For any finite time series there is a sufficiently low 
frequency such that the magnetization switches during 
every half-period of the field. For very low frequencies 
the magnetization switches before the field reaches its ex- 
treme value during every half-period. In this frequency 
regime, m(t) switches sign abruptly relative to the length 
of a field period (see Fig. ||(d)). The switching time t s 
may be found by solving m(t s ) = 0. Thus the equation 



In 2 = $(i„ 



(6.1) 



As Figs. ||(d) and ||(d) show, for sufficiently low fre- 
quencies m(t) switches sign early in the period, i.e. at 
a value |i?(t)| < Hq. This allows the approximation 
H{t) « Hocut to be used for the field, henceforth re- 
ferred to as "the linear approximation in the field." Us- 
ing this approximation and the substitutions x — H^Lot 
and x' = HQU>t n allows Eq. (|5.3|) to be rewritten as 



2 r H s 



B(T)tt 2 v 
B{T)n 2 v 2 ( 



x' 3 (x 2 -x' 2 ) 2 exp 



dx' 



x I x' 3 exp 



dx 1 



-2x 2 I x' 5 exp 



socn 



dx' 



(6.2) 



'7 

x exp 



socn 



dx' 



where H s — H^ujt s . Together with Eq. (6.1) this yields 



2H*E b (T)T 



-4,- 



H. 



+El(T)T 



-8,- 



where each of the integrals in Eq. ( p.2[) has been rewritten 
as an incomplete gamma function |89[ using the expres- 
sion 




(6.3) 



u e 



- a ' u du 



,n+l 



r 



(6.4) 



For small w the hysteresis loops are practically square, so 
the scaled loop area in the low-frequency (LF) limit can 
be expressed as 



LF 



iH 



Ho 



(6.5) 



The switchi ng fi eld H s (lu) is obtained by numerical solu- 
tion of Eq. (6.3), and the result for (^4) LF /4-ffo is shown 
as the dotted curve in Fig. |^. There is good agreement 
between the linear approximation calculation, the numer- 
ical integration calculation (NI) , and the MC data for low 
frequencies. The slight overestimate of the loop area by 
the linear approximation at low frequencies is due to a 
systematic error in the way that the loop area is calcu- 
lated from the value of H s (uj). The disagreement near 
the maximum in the loop area is due to the breakdown 
of the linear approximation as the magnetization begins 
to switch sign only at fields close to Hq. 
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The dash-dotted curve in Fig. ^ is the theoretical low- 
frequency prediction for L = 64 in the SD regime. This 
curve corresponds to the solid curve in Fig. 11(b) of 
Ref . [p(| , appropriately re-scaled so that it may be shown 
together with the MD results. This re-scaling consists of 
using Hq = 0.3J in the quantity A/4Hq and (t(Ho)) for 
the MD regime in calculating the scaled frequency 1 / R. 
The MC data point at the very lowest frequency in Fig. ^ 
(i? = 5000) agrees with the SD, rather than the MD, the- 
oretical prediction, even though the simulation was per- 
formed with H = 0.3 J. This crossover from MD to SD 
hysteretic behavior is a completely frequency-dependent 
effect. As the frequency of the field decreases, the value 
of the switching field, H s (u>) decreases as well. For suf- 
ficiently low frequencies, H s (oj) < Hbsp{L,T), the sys- 
tem undergoes SD decay before the field becomes large 
enough to produce MD decay. In fact, the intersection of 
the SD and MD theoretical predictions occurs at a value 
of A which corresponds to a square loop with a switching 
field of H s (uj) m H DS p(L,T). While m(i) and A do not 
depend on system size in the MD regime, A in the SD 
regime, and hence the location of the crossover, depends 
on L through the L~ d dependence of the lifetime in the 
SD regime (see Eqs. (4.15) and (7.13) of Ref. |oJ). 

One can obtain an approximate analytic solution for 
(A)t F by taking the first three terms in the asymptotic 
expansion fl89| 



T(a, x) 



1 (a-l)(a-2) 



(6.6) 



for large x. Straightforward, but lengthy, a lgeb ra gives 
the following asymptotic expansion for Eq. (6.3) 



So(T) 



exp 



where we define 



D = 



socn 



In 2 



2B(T)Sl 2 v 2 E*(T) 



(6.7) 



(6.8) 



With B(T) = 0.02048J- 3 MCSS _1 and the values found 
in Table § this gives D = 17.905 J^MCSS" 1 . Ignoring 
the nonexponential prefactor in Eq. ( |6.7|) , solving fo r the 
switching field H s (lo), and substituting it into Eq. (6.5) 
results in the asymptotic, logarithmic frequency depen- 
dence for the loop area, 



(A) 



LF 



-3o(T) [-HDH io)Y 



(6.9) 



As in the SD case j5^,[52| , from a log-log plot of the loop 
area versus frequency one can extract effective exponents 
from the data over nearly two decades in 1/R. How- 
ever, these effective exponents depend strongly on the 
frequency range in which the fit is performed. Very small 



effective exponents that appear consistent with this pic- 
ture were recently reported by Suen et al. for ultrathin 
Co films on Cu(001) @. 

As for the SD regime, we stress that the asymptotic 
low-frequency behavior described by Eq. (6.9) would only 
be seen for extremely low frequencies, where the dotted 
and dashed curves in Fig. |] come together. Therefore, 
due to the dependence of 7?dsp on L the crossover to SD 
droplet behavior will occur at frequencies much larger 
than those for which this asymptotic expression is valid 
even for large L. For a system size of L = 10 9 the ratio of 
the full numerical solution to the asymptotic expression 
is approximately two near the crossover frequency, which 
for that system size would be near 1/R ss 6 x 10~ 9 . 

Relatively few studies have considered B, the correla- 
tion of the magnetization with the external field. Our 
theoretical derivations of B are analogous to those for 
A. The probability densities of B in the MD regime 
are shown in Fig. [l(]. Figure |l(](a) shows the probabil- 
ity densities for several values of R between 2 and 200. 
The source of the roughness of the distributions at the 
lowest frequencies is the same as in Fig. |?](a). Figure 
|l^(b) shows the probability densities for a frequency near 
the transition for L = 64, 90, and 128. The system size 
dependence of the peak position and the width of the 
distributions is similar to that for the loop-area distribu- 
tions. The means of the distributions in Fig. ^Q(a) are 



shown in Fig. |ll] along with the results for L = 90 and 
128. T he so lid curve is obtained by numerical integration 
of Eq. Jl.4j ) similar to the procedure for A. The agree- 
ment between the MC data and the theoretical predic- 
tion is very good for low and high frequencies, but poor 
for intermediate frequencies, where the MC data take on 
negative values. This disagreement is probably due to 
the procedure for calculating m{t) during a half-period 
in which shrinking droplets are present in the system (see 
Sec.0). 

As for the SD regime (|(]], the physical significance of 
the integrals A and B can be clarified by comparison 
with linear response theory. One can easily find that 
A/(ttHq) and 2B/Hq correspond to the dissipative and 
reactive parts of the complex linear response function, 
respectively. It is therefore natural to combine A and B 
into an analogous nonlinear response function, 



A(ff ,7» = 



1 

IP- 



2B 



-A 



(6.10) 



The frequency dependence of the norm of this response 
function, \X\, is shown in Fig. [l^. The maximum in A 
(Fig. ||) and the sign change in B (Fig. [ll]), which occur 
close together in frequency, are characteristic behaviors of 
dissipative and reactive parts of a response function near 
resonance. The behavior of (A) and (B) in this frequency 
range is very similar to that observed in the SD regime 
|50fl, where we associated it with stochastic resonance 
|90|. Whether or not this name is appropriate for the 
low-frequency synchronization of m(t) with H(t) in the 
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MD regime, is probably a matter of taste. Although the 
overall fluctuations in m(t) are small, the switching is 
entirely driven by random nucleation on a microscopic 
scale. 



VII. DYNAMIC PHASE TRANSITION 

Although nonequilibrium phase transitions have been 
studied for over two decades, the understanding of their 
universality and scaling properties remains much weaker 
than for equilibrium critical phenomena. Finite-size scal- 
ing is routinely used to determine the location and critical 
exponents of equilibrium, thermodynamic phase transi- 
tions from numerical results for finite systems . A 
true phase transition occurs only in the thermodynamic 
limit, i.e. as the system size approaches infinity while the 
energy and particle densities are kept constant. For sim- 
ulations on finite systems, quantities such as the suscep- 
tibility and specific heat only display a bounded peak, 
rather than a divergence, as a function of temperature 
or field. However, for simulations performed on larger 
and larger lattices, thermodynamic quantities approach 
the infinite-size limit. It is this approach to the thermo- 
dynamic limit that is utilized to obtain estimates of the 
critical exponents from various system sizes. 

Second-order phase transitions in equilibrium systems 
are characterized by a set of critical exponents, each 
of which describes the behavior of a different quantity 
at the critical point. Three of these quantities and 
their associated exponents are [ p2[ : the order parameter 
(M - \if for i < 0), the susceptibility ( X ~ |i|~ 7 ), and 
the correlation length (£ ~ where ~ denotes "the 

singular behavior of." The field conjugate to the order 
parameter is understood to be zero, and i= (T — T c )/T c 
is the reduced temperature. Finite-size scaling theory al- 
lows one to estimate the critical exponents by measuring 
the system size dependence of various quantities. Com- 
bining the expression for the correlation length exponent 
with the finite-size scaling assumption, £(T C (L)) ~ L, 
gives 



\T C (L) - T c \ cx L 



-l/u 



(7.1) 



where T C (L) can be defined as the location of the peak in 
the susceptibility for a given L W^J&lti. When combined 



with the definitions for the critical exponents, Eq. (7.1) 
gives 



and 



^poak a L7 /„ 



(\M\ n ) L cx L" n{p/v) 



(7.2) 



(7.3) 



where x? 6 is the maximum value of the susceptibility 
for a given L, and (|M| n ) L is the nth moment of the 
norm of the order-parameter at T c . 



The period-averaged magnetization Q has been pro- 
posed as a "dynamic order parameter" for systems ex- 
hibiting hysteresis |35|,|7]-|2| 69 1 . Those studies of the 
Ising model have suggested the existence of a dynamic 
phase transition between an ordered dynamic phase with 
(\Q\) > and a disordered dynamic phase with (|Q|) w 0. 
Figure [l3| shows the probability densities of Q in the MD 
regime for L = 64, 90, and 128. For each system size, 
as the frequency of H(t) decreases, the probability den- 
sities for Q change from bimodal distributions with the 
two peaks each centered around a nonzero value of Q, to 
unimodal distributions with a peak around Q = 0. Each 
of the Q distributions shown in the three-dimensional 
Fig. [l3| is a histogram of Q time-series values (see Fig. ||) 
at a particular frequency. These distributions suggest the 
presence of a second-order dynamic phase transition. In 
particular, the frequency dependence for these probabil- 
ity densities is strikingly similar to the dependence on 
inverse temperature for probability densities of the equi- 
librium magnetization in the zero-field Ising model. We 
therefore identify the norm of the period-averaged mag- 
netization, \Q\, as the order parameter of the dynamic 
phase transition, and we apply finite-size scaling theory 
in analogy to the scaling theories used to quantify second- 
order phase transitions in equilibrium systems. Figure [l4| 
shows the average norm of the period-averaged magne- 
tization, for the same system sizes as in Fig. |lj. 
This figure clearly suggests a DPT as the average order 
parameter (\Q\) changes from a value near zero to a non- 
zero value. Rather than a sudden change in the order 
parameter, the transition region is "smeared" out due to 
finite-size effects. The mean and standard deviation for 
\Q\ are system size dependent as well. We quantify these 
finite-size effects below. 

At a second-order phase transition there is a diver- 
gence in the susceptibility. For equilibrium systems, the 
fluctuation-dissipation theorem relates the susceptibility 
to fluctuations in the order parameter. For the present 
system, it is not obvious what the field conjugate to Q 
might be. Therefore, we cannot measure the susceptibil- 
ity directly. However, we can calculate the variance in 
\Q\ as a function of frequency and study its system size 
dependence. We define X as 



X = L 2 Var(|Q|) = L 2 (Q 2 ) - (|Q|) 



(7.4) 



If the system were to obey a fluctuation-dissipation rela- 
tion, X would be proportional to the susceptibility, and 
both would scale with L in the same manner. Figure 
[l5| shows X vs 1/R for all three system sizes. For all 
three values of L, X displays a prominent peak near the 
transition frequency, which increases in height with in- 
creasing L, while no finite-size effect is seen at lower and 
higher frequencies. This finite-size effect in X implies 
the existence of a divergent length associated with the 
order-parameter correlation function near the dynamic 
transition. The observation that P(|Q|) displays no peak 
near \Q\ =0 in the ordered dynamic phase is additional 
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evidence for the second-order (as opposed to first-order) 
nature of this transition |}3) . 

The cumulant intersection method (7||)l|] is useful 
for determining the location of a second-order transition 
when the critical exponents are not known. In order to 
estimate the location of the transition wc define the "dy- 
namic" fourth-order cumulant ratio 



U L = 1- 



3<IQI 2 )i 



(7.5) 



where (|Q| n ) = /„ \Q\ n P{\Q\)d\Q\. Figure |l6 shows U L 
vsl/R for the same system sizes shown in Fig. 15]. Above 
the transition frequency, in the > ordered dy- 

namic phase, Ul approaches 2/3, corresponding to two 
narrow peaks centered at ± Below the transition 

frequency, in the ~ disordered dynamic phase, 

Ul approaches 0, corresponding to a Gaussian centered 
at zero. At the transition, the cumulant should have a 
non-trivial fixed value, U*. Therefore, the location of 
the cumulant intersection gives an estimate of the transi- 
tion frequency without foreknowledge of the critical expo- 
nents. Due to the large spacing of our data and possible 
correction-to-scaling effects, we cannot identify a unique 
intersection point. We estimate the location of the inter- 
section by the crossing for the two largest system sizes 
near R' 1 m 0.2910 (i? cr ps 3.436) with U L = U* m 0.61. 
This is close to an extremely precise transfer matrix cal- 
culation of U* = 0.6106901(5) @ as well as MC esti- 
mates p5| for the two-dimensional Ising model. How- 
ever, as recently pointed out by Luijten et al. |96|, the 
value of the cumulant intersection should not be taken 
too seriously unless a sufficient range of system sizes is 
available. 

With our estimate for the transition frequency, we can 
now approximate the critical exponent s /3, 7, an d v char- 
acterizing the transition by using Eqs. (7.1)-(f7~5|), replac- 



ing T c by R c , X *l by , and (\M\ n ) L by (\Q\") L . 

To extract exponent estimates using these relations, we 
use a method sometimes referred to as "phenomenologi- 
cal renormalization" [glj of the MC data. This method 
consists of estimating an exponent by using two system 
sizes, bL and L. The following example is a derivation of 
an exponent estimate for [3 /v. 

(&£)-*»(0/«O 



(\Q\ n ) 
(\Q\ n ) 

which yields 

~(\Q\ n ) bL 



bL 



hi 



(\Q\ T 



In b 







+ 0(1/ Info) . (7.6) 



Similar relations can be found for the other exponent 
ratios, 




lnfo = ~+0(l/ln&) 



lnfo = - + <3(l/lnfo) 
v 



(7.7) 
(7.8) 



In the large-system limit these exponent estimates will be 
linear when plotted vs (Info) -1 . Then, one can extrapo- 
late to the infinite size limit by performing a linear fit of 
the data to find the intercept at (Info) -1 =0. Simulations 
with larger system sizes would be computationally pro- 
hibitive, and smaller system sizes would no longer be in 
the MD regime. With data for only three system sizes, 
the exponent estimates obtained using the two largest 
system sizes are easily shown to be identical to those ob- 
tained using the extrapolation procedure above. We cal- 
culate two sets of estimates for (3/v, one using the scaling 
relation for the second moments of the order-parameter 
distribution (n = 2) and the other using (n = 4), ob- 
taining {P/v) n =2 ~ 0.111 and {(3/v) n =A ~ 0.113. Our 
estimates for the other exponents are j/v f=a 1.84 and 
v w 1.1. Also, we obtained an independent estimate for 
the exponent v by measuring the finite-size effects in the 
location of the high-frequency zero-crossing in (B). The 
estimate obtained is v « 1.09, in good agreement with 
that obtained from the location of X poak . Our results are 
close to the two-dimensional Ising values for the analo- 
gous exponent ratios (/3/^ = 1/8 = 0.125, 7/f = 7/4 = 1.75, 
v= 1). Given the accuracy of our data however, our ex- 
ponent estimates are also not inconsistent with the uni- 
versality class of two-dimensional, random percolation 
{P/v = 5/48 w 0.104, j/u = 43/24 « 1.79, !/ = 4/3« 1.33). 
Combining the exponent estimates we find 



2(PM + (iM « 2.06 : 



(7.9) 



This relation between the measured exponent ratios indi- 
cates the consistency of our scaling procedures and thus 
strengthens our belief that the dynamic transition is a 
genuine, continuous phase transition. If the divergent 
length is indeed the correlation length that descri bes the 
order-parameter correlation function, then Eq. (7.9) is 
a hyperscaling relation. The DPT critical point then 
should be a nontrivial fixed point in the renormaliza- 
tion group sense. Based on the evidence presented in the 
following paragraph, we believe this is the case. 

One should also consider the possibility that hyper- 
scaling is violated and the DPT represents a mean-field 
critical point. The divergent length would then be the 
"thermal length" |7^,|9(|, whose divergence is governed 
by the exponent (2/3 + 7)/d. This wo uld then be the 
exponent we have called 'V," and Eq. (7.9) would hold 
exactly as a tautology. However, our estimates for (3 and 
7 are far from those of a mean- field 4 model (/3 = 1/2 
and 7 = 1). Likewise, our estimated cumulant crossing, 
U* ~ 0.61, is far from the expected mean-field value, 
U* = 1 - r 4 (l/4)/(247r 2 ) w 0.27 @f7|. Furthermore, 
if the phase transition were to be mean-field in this two- 
dimensional system, it should have to be induced by some 
effective long-range interaction, which then should have 
the same effect in one dimension. However, exploratory 
MC simulations indicate that the one-dimensional Ising 
model in an oscillating field does not have a dynamic or- 
dered phase ]9q| . The evidence summarized in this para- 
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graph makes it extremely unlikely that hyperscaling is 
violated by the DPT. 

The consistency of the estimates of v from the posi- 
tions of X pcak and the high-frequency zero-crossing of 
(B) indicates that this zero occurs at the DPT. The two 
zeros of (B) are clearly separated in frequency, and our 
finite-size scaling results indicate that they remain so as 
L — > oo. The low- frequency zero is associated with the 
maximum in (A) . These observations enable us to answer 
the question recently raised by Acharyya |8^] of whether 
the DPT corresponds to the maximum in (A). It does 
not. 

To further illustrate the nature of the dynamic phase 
transition, the finite-size effects in the distributions for 
the norm of the order-parameter, |Q|, are shown in 
Fig. [I?]. These probability densities for \Q\ are the 
symmetrized versions of selected distributions shown in 
Fig. [l3l Distributions in the ordered dynamic phase re- 
gion, i.e. above the transition frequency, should move to- 
ward a constant, nonzero value of \Q\ and become nar- 
rower with increasing L. This is seen in Fig. |l^(a). For 
this frequency, the distribution of Q for L = 90 is highly 
asymmetric about zero and for L = 128 the distribution is 
unimodal. This gradual loss of symmetry with increasing 
L is due to the finite length of the simulation time series, 
but it does not adversely affect our ability to analyze 
P(|Q|). The distributions in Fig. |l7](b) are in the disor- 
dered dynamic phase region, i.e. at a frequency slightly 
below the transition frequency. Due to finite-size effects 
however, the distributions for L = 64 and 90 appear to be 
centered about nonzero values of \Q\. The distributions 
in Fig. [l7](c) are near the transition and should scale with 
system size L. We assume that the mean of the order 
parameter scales with L and define the scaling variable 
Q = L/ > '"\Q\. Hence, the scaled probability density for 
\Q\ is given by 

P L (Q) = L-?/»P(\Q\) , (7.10) 

where the prefactor L~^l v ensures conservation of prob- 
ability. Figure |l8| shows this scaled probability density. 
The peak positions scale fairly well, the peak heights less 
so. This could be due to the following reasons. The fre- 
quency might be sufficiently far from the transition that 
single-parameter scaling is not adequate, and there might 
be corrections to the finite-size scaling that are large for 
these relatively small system sizes. Also, the lack of scal- 
ing for the peak heights could be due to the asymmetry 
in P(Q) near the transition. 

The results in this section clearly show that the statis- 
tical properties of the order parameter Q exhibit finite- 
size scaling, and that scaling techniques developed for 
estimating the critical exponents for second-order phase 
transitions in equilibrium systems apparently can be suc- 
cessfully applied to estimate the exponents associated 
with the dynamic phase transition. While these scal- 
ing relations are concerned with \Q\, it is worth men- 
tioning that one may also measure the fluctuations in 



the other two quantities measured, A and B. Figure |l^ 
shows the fluctuations for A and B, defined in analogy 
to the order-parameter fluctuation X. The fluctuations 
in A seem to show slight finite-size effects, as the peak 
positions appear to be approaching i?" 1 with increasing 
L. One might speculate that this could indicate that A is 
coupled to energy fluctuations which are logarithmically 
divergent as they are for the two-dimensional Ising model 
in equilibrium. The fluctuations in B show no significant 
finite-size effects. 



VIII. DISCUSSION 

The mechanism by which a metastable phase decays 
depends sensitively on the system size, the temperature, 
and the strength of the applied field. For large sys- 
tems and moderately strong fields, the decay proceeds 
through the nucleation and growth of many droplets of 
overturned spins in different parts of the system. This 
regime has been termed the multi-droplet (MD) regime. 
In this regime the magnetization response in static field is 
described by the KJMA approximation (Avrami's law), 
which assumes the presence of many noninteracting, 
overlapping droplets. Theoretical predictions by a gener- 
alization of the KJMA approximation, in which a time- 
dependent nucleation rate and droplet interface velocity 
are used, agree well with simulations for quantities like 
the average hysteresis-loop area and correlation, espe- 
cially at low driving frequencies. The time dependence 
is included in the theory by replacing the constant field 
H by H(t) = — Hq s'm(uit). This central idea provides 
the analytic framework for our theoretical descriptions of 
quantities measured from the MC simulations in the MD 
regime. The theoretical calculations and the MC data 
agree very well, especially considering that only one ad- 
justable parameter is needed, which was measured from 
a particular hysteresis simulation (R = 200). All of the 
other constants used are either known from droplet the- 
ory or were measured for MC simulations of field reversal 
in kinetic Ising models. To the best of our knowledge, the 
work reported here is the first which explicitly considers 
hysteresis for the Ising model in the MD regime. 

We compute the power spectral densities from the sim- 
ulated time series and qualitatively explain various fea- 
tures of the spectra in the full frequency range from the 
lowest observable frequencies to the rapid fluctuations 
due to thermal noise. For low field frequencies, the sys- 
tem is in the disordered dynamic phase, and the time 
series contain no large fluctuations. Consequently, the 
PSDs are flat at frequencies below the fundamental peak 
at the frequency of the field. The significant power den- 
sity in the low-frequency portion of the PSDs corresponds 
to the long-time behavior in the filtered time series for 
Q. For high field frequencies, the system is in the ordered 
dynamic phase, and the time series display long-time be- 
havior as the system switches between thermodynamic 
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phases. This corresponds to large power density in the 
low-frequency portion of the PSDs. Near the dynamic 
phase transition the PSDs exhibit similar behavior in the 
low-frequency part of the spectrum. Due to insufficiently 
long time series we are unable to resolve any quantitative 
difference between these PSDs and those in the ordered 
dynamic phase region. 

We also calculate the hysteresis-loop area A and the 
correlation B for a wide range of frequencies. Because 
of its role as a measure of the energy dissipation in the 
system, A is a quantity of particular experimental signif- 
icance. For all frequencies, the loop-area and correlation 
distributions are unimodal due to the almost determinis- 
tic magnetization response in the MD regime. Our the- 
oretical predictions for the frequency dependence of (A) 
and (B) use the time-dependent extension of Avrami's 
law to calculate m(t), from which the loop area and cor- 
relation are calculated explicitly. The assumption is that 
the m(t) values calculated for a single period accurately 
describe the average values of A and B over a long sim- 
ulation time series. This assumption should be expected 
to break down most significantly for frequencies near the 
dynamic phase transition, where the fluctuations in the 
magnetization response are largest. This is clearly seen 
in Figs. H and [ll], where the least satisfactory agree- 
ment between the theory and the MC data occurs for 
frequencies near the dynamic phase transition. For the 
low-frequency regime we obtain an analytic expression 
for (A) . Our theoretical calculation agrees well with our 
MC results and predicts an extremely slow crossover to 
a logarithmic dependence of the loop area on H uj. The 
switching dynamics is dominated by nucleation and indi- 
cates no overall power-law dependence for the loop area 
on field amplitude and/or frequency, in contrast to what 
has been claimed in other simulational and experimen- 
tal studies. We emphasize that numerical analysis of data 
generated by our analytic solution, even over two or three 
frequency decades, could easily lead to the conclusion 
that the data were taken from a power law. Our simula- 
tions reveal that for frequencies far higher than those at 
which the asymptotic logarithmic dependence would be 
observable, a system-size dependent crossover from MD 
to SD behavior occurs. This novel frequency dependence 
for A is a consequence of the field dependence of the SD 
and MD decay mechanisms. As the frequency of the field 
becomes sufficiently small, the system is subject to fields 
smaller than Hdsp (T, L) for a sufficiently long time so 
that SD decay usually occurs before the field becomes 
large enough for MD decay to happen. 

The period-averaged magnetization Q has been pro- 
posed as an order parameter associated with the dy- 
namic phase transition (DPT) in kinetic Ising models. 
The DPT is a nonequilibrium phase transition which oc- 
curs due to an explicit time-dependence in the Hamil- 
tonian, rather than the dynamical rules governing the 
system. The probability densities that we obtain for Q 
clearly show that the system changes from an ordered dy- 
namic phase with nonzero to a disordered dynamic 



phase with (|Q|) ~ as the field moves from high to 
low frequencies. To distinguish this frequency-dependent 
change in (\Q\) as a true second-order phase transition 
rather than merely a simple bifurcation, we measure the 
finite-size effects at the DPT and apply finite-size scal- 
ing (FSS) techniques analogous to those used to mea- 
sure the critical exponents which characterize equilib- 
rium second-order transitions. The measured exponents 
(fi/v ~ 0.11, "f/v m 1.84, and v w 1.1) are close to both 
the two-dimensional Ising and random percolation val- 
ues, and they represent strong evidence that hyperscal- 
ing is obeyed. Our success in applying FSS techniques 
borrowed from the theory of equilibrium second-order 
phase transitions to this nonstationary nonequilibrium 
problem suggests the possibility of mapping other suit- 
ably defined quantities for this system to thermodynamic 
entities, such as the field conjugate to the order parame- 
ter, the specific heat, and the correlation length. Such a 
nonequilibrium thermodynamic theory for steady states 
as been attempted by Paniconi and Oono 

In future work we plan to analyze longer simulations 
on larger system sizes to more accurately determine the 
values of the exponents and the location of the DPT. This 
includes measuring the possible finite-size effects in the 
probability distributions of the energy |65[ , which might 
be related to the finite-size effects seen in the fluctua- 
tions in A. If a fluctuation-dissipation theorem could be 
shown for this system, the order-parameter and energy 
fluctuations could be related to a nonequilibrium suscep- 
tibility and specific heat respectively. Another important 
question left to future study is if and to what extent the 
exponents depend on the temperature and field ampli- 
tude. While the critical frequency will almost certainly 
depend on amplitude and temperature, the critical expo- 
nents would most likely not if the DPT indeed represents 
a new "dynamic universality class." 

Finally, we note that the quantities that we have ana- 
lyzed numerically could all be measured in experiments 
on hysteresis in a variety of systems and analyzed by 
methods essentially identical to our analysis of the MC 
data. 
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TABLE I. Parameters and constants used in this work. The values of the parameters Ho, L, and T have been selected such 
that switching occurs via the multi-droplet mechanism. The constants So (T) and K are calculated from droplet theory Jt4| [rjj] 
for two-dimensional Ising systems. The constants (r) and r are measured from field-reversal MC simulations with the Glauber 
dynamic (using the parameters listed in the left column). The constants fi(T) j77| and v (T) frjj], where the droplet interface 
velocity is Vo = v(T)\H\, have been measured in other work. The value for J/dsp is taken from Fig. 11 of Ref. fed]. For L = 90 
and 128, the relative standard deviation is r — 0.072 and 0.053 respectively. Except for -Hdsp, which decreases slowly with 
increasing L, all other values are the same for L — 90 and L = 128. 



Parameters 


Constants (theory) 


Constants (simulation) | 


Ho 


0.3J 


So(T) 


0.506192 J 


Oa(T) 


3.15255 


L 


64 


K 


3 (exact) 


u(T) 


(0.465 ± 0.014) J^MCSS" 1 


T 


0.8T C 




(r) 


74.5977 MCSS 






Hdsp (L = 64) 


(0.11 ±0.005) J 






r(L = 64) 


0.105 



0.35 
0.3 

0.25 
0.2 

0.75 

a j 

0.05 
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FIG. 1. Location of MD simulations in the (Ho,L) plane. The solid curve represents the dynamic spinodal (DSP), 
Hdsp ~ (In L)~ 1 ^ d ~ 1 K This theoretical curve is an asymptotic result obtained by setting Rq ~ L. The filled circles denote the 
system sizes (L = 64, 90, 128) and field amplitude (J/o = 0.3J) used here to study hysteresis in the MD regime. The open circle 
denotes the system size and field amplitude used to study hysteresis in the SD regime in Ref. [|o| . The dotted curve represents 
a theoretical result for the "thermodynamic spinodal field" (THSP) p7| , ^8| which separates the SD region from the coexistence 
regime (CE). It is obtained by setting R c w L. 
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(b) R = 3.436 




(d) R = 200 



10000 20000 30000 40000 

FIG. 2. Short segments of the magnetization time series 

m{i) (solid curve) and the external field H(t) (dashed curve) 

vs timet forT = 0.8T c , d = 2, L = 64, and .Ho = 0.3J. The total 

length of the time series is approximately 16.9 x 10 6 MCSS. 

For these parameter values the average lifetime in static field 

is (t(Ho)) ~ 75 MCSS. The time series are shown for the 

scaled field periods (a) R=3, (b) R = 3.436 « R CI , (c) R=7, 

and (d) i? = 200. 
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(b) R = J.4J6 






J 






(e) R=200 




0.5 



















^- ~— - 




— 


"- 




-0.5 










-1 















10 


20 

n^ritulv 


30 



FIG. 3. Short segments of the "filtered time series" of 
period-averaged magnetization values, Q, vs number of field 
periods. The parameters used are the same as in Fig. |^. (a) 
R = 3, (b) R = 3.436 « R CI , (c) R = 4, (d) R = 7, and (e) 
i?=200. 
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FIG. 4. (a) Power spectral densities (PSDs) for L = 64. 
Spectra are shown for three different frequencies of the field, 
and are plotted with an arbitrary offset for clarity. In addition 
to a change in the amount of smoothing, the right-hand sec- 
tion of each spectrum contains only one out of every 250 data 
points to facilitate plotting. The magnetization is sampled 
every 1.0 MCSS, so the Nyquist frequency is (On/27t) =0.5 
MCSS -1 . The lowest frequency that can be resolved is 
2.38 x 10~ 7 MCSS -1 . The dashed line with slope -2 is a 
guide to the eye. (b) Same spectra without the offset to il- 
lustrate how all three PSDs fall near the thermal noise back- 
ground at high frequencies. Note that the spectra for R = 3 
and R = 3.436 « R CI cross at low frequencies. 




FIG. 5. Schematic of growing and shrinking droplets in 
the MD regime for a sufficiently high frequency such that 
m(t) does not completely switch during a period. In (b) and 
(c) the dark regions represent the stable phase, and the light 
regions represent the metastable phase. The arrows indicate 
the growth direction of the droplet interfaces. 
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FIG. 6. Representative hysteresis loops obtained from 
Monte Carlo simulation data with system size L = 64 for (a) 
R=3, (b) R = 3.436 » R CI , (c) R=7, and (d) i?=200. Each 
panel shows loops for the same time intervals shown in the 
corresponding time-series data in Fig. ^. 
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FIG. 7. Probability densities of the hysteresis-loop area, 
A = — <f> rn(H) dH. The loop area has been normalized by 
the maximum possible loop area, 4Ho. (a) L = 64. The values 
of the scaled period shown are R = 2, 3, 3.436 w _R cr , 3.5, 
3.75, 3.9, 4, 5, 6, 7, 8.4, 12, 25, 50, 80, 140, and 200, so 
the R axis is not linear. The two R values explicitly marked 
indicate the directions of increasing and decreasing frequency, 
(b) Distributions near the critical frequency shown for L = 64, 
90, and 128. The distributions are narrower away from R CI 
but look qualitatively similar. 
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FIG. 8. Mean and standard deviation of the loop-area 
distributions vs the scaled frequency 1/R. The data points 
are the means of the distributions shown in Fig. Q for L = 64 
along with the means of the corresponding distributions for 
L = 90 and 128. The vertical bars are not error bars, but 
give the standard deviations of those distributions. The solid 
curve comes from numerical integration of Eq. ( [l.3] ) using the 
values of m(t) from Sec. The single free parameter, B(T), 
is adjusted so the theoretical prediction agrees with the data 
point at 1/J?=0.005 as described in the text. 
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b S ,oU/R) 

FIG. 9. Log-log plot for the mean of the loop-area distri- 
butions vs the scaled frequency 1/R. The vertical bars are 
not error bars, but give the standard deviations of those dis- 
tributions. The data points are the same as those used for 
L — 64 in Fig. |^ for the lowest frequencies. The solid curve 
(NI) is the same as the solid curve in Fig. ^. (Due to nu- 
merical difficulties, this calculation was not extended to lower 
frequencies than those shown.) The dotted curve re sult s from 
numerical solution of the linear approximation, Eq. ( |6.3|) . The 
dash-dotted curve is the prediction for the loop area in the SD 
regime for L — 64, scaled so that it may be plotted along with 
the MD results. The dashed curve is the asymptotic, logarith- 
mic frequency dependence for the MD loop area, Eq. (6.9). 
This asymptotic result approaches the full solution only for 
frequencies that are lower than the crossover to the SD solu- 
tion, even for very large L. The arrow indicates the area of a 
hysteresis loop with H s « _Hdsp(£ = 64). 
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FIG. 10. Probability densities of the correlation, 

B = (w/2tt) §m{t) H(t) dt. fa) L = 64. The values of R 
shown are the same as in Fig. W. (b) Distributions near the 
critical frequency shown for L = 64, 90, and 128. The lines are 
guides to the eye. The distributions look qualitatively similar, 
even away from R CI . 
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FIG. 11. Mean and standard deviation of the correlation 
distributions vs the scaled frequency 1/R. The data points 
are the means of the distributions shown in Fig. 0(a) for 
L = 64, along with the corresponding results for L — 90 and 
128. The vertical bars are not error bars, but give the stan- 
dard deviations of the distributions. The solid curve comes 
from a theoretical calculation analogous to that in Fig. ^ 
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FIG. 12. Mean and standard deviation of the norm of the 
response function, \X\, vs the scaled frequency 1/R. The data 
points are obtained using the loop-area and correlation data 
in Fig. |^ and Fig. [H]. The solid curve is obtained using the 
theoretical values for (^4) and (B) in those figures. 
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FIG. 14. Mean of the norm of the period-averaged mag- 
netization vs the scaled frequency 1/R. The finite-size effects 
are clearly seen for frequencies in the neighborhood of the dy- 
namic phase transition. The arrow indicates the approximate 
value of the critical frequency 1/R CI . Lines connecting the 
data points are guides to the eye. 
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FIG. 13. Probability densities of the period-averaged 
magnetization, Q = (lu/2tt) <j> m(t) dt. (a) L = 64. The values 
of R shown are the same as in Fig. @. (b) L = 90. The values 
of the scaled period shown are R = 3, 3.436 w _R cr , 3.5, 3.75, 
3.835, 4, and 5. (c) L — 128. The values of the scaled period 
shown are R = 3, 3.436 w R C r, 3.5, 3.612, 3.693, 3.721, 3.75, 
4, 5, 7, and 200. The unimodal distributions for the smallest 
R value in (a) and (c) are due to the finite simulation time. 
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FIG. 15. The order-parameter fluctuation 

X — L 2 Var(|Q|) vs the scaled frequency 1/R. The "disor- 
dered dynamic phase," ((|Q|) ~0), lies on the low- frequency 
side of the peaks. The "ordered dynamic phase," ((|Q|) > 0), 
lies on the high-frequency side. The statistical error bars are 
estimated by partitioning the data into ten blocks. Error bars 
smaller than the symbol sizes are not shown. The arrow in- 
dicates the approximate value of the critical frequency 1 /R CI . 
Lines connecting the data points are guides to the eye. 
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FIG. 16. Fourth-order cumulant ratio Ul vs scaled fre- 
quency 1/R, for L = 64, 90, and 128. We use the same sym- 
bols as in Fig. [l5| The horizontal line marks (7l = 2/3. Lines 
connecting the data points are guides to the eye. Inset: area 
close to the cumulant crossing at 1/R C1 - ~ 0.2910. 
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FIG. 17. Probability distributions for the norm of the pe- 
riod-averaged magnetization \Q\ for a frequency (a) above the 
transition, 1/R — 0.333, (b) below the transition, 1/7? = 0.25, 
and (c) near the transition 1/7? = 0.2910. Lines connecting 
the data points are guides to the eye. 
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FIG. 18. Scaled probability densities of \Q\ for three sys- 
tem sizes. The same symbols are used as in Fig. The 
scaling function is L~®'" P(\Q\) vs L^'"\Q\, and the value of 
the scaling exponent used is (f5/v) n =2 ~ 0.11. The scaled fre- 
quency of the field is 1/R = 0.2910 « l/R C i- Lines connecting 
the data points are guides to the eye. 
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